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3.6 Solving Large Linear Systems 

Finite elements and finite differences produce large linear systems KU = F. The 
matrices K are extremely sparse. They have only a small number of nonzero entries in 
a typical row. In "physical space" those nonzeros are clustered tightly together — they 
come from neighboring nodes and meshpoints. But we cannot number nodes in a 
plane in any way that keeps neighbors close together! So in 2-dimensional problems, 
and even more in 3-dimensional problems, we meet three questions right away: 

1. How best to number the nodes 

2. How to use the sparseness of K (when nonzeros can be widely separated) 

3. Whether to choose direct elimination or an iterative method. 

That last point will split this section into two parts — elimination methods in 2D 
(where node order is important) and iterative methods in 3D (where preconditioning 
is crucial). 

To fix ideas, we will create the n equations KU = F from Laplace's difference 
equation in an interval, a square, and a cube. With unknowns in each direction, 
K has order n = N or N"^ or N^. There are 3 or 5 or 7 nonzeros in a typical row of 
the matrix. Second differences in ID, 2D, and 3D are shown in Figure 3.17. 



Tridiagonal K 



-1 # Block 

Tridiagonal 



4 



K 




by N"^ 

Nhy N -1 4 by A^^ 



Figure 3.17: 3, 5, 7 point difference molecules for —Uxx, —Uxx — Uyy, —Uxx — Uyy — u 



Along a typical row of the matrix, the entries add to zero. In two dimensions this 
is 4 — 1 — 1 — 1 — 1 = 0. This "zero sum" remains true for finite elements (the element 
shapes decide the exact numerical entries). It reflects the fact that u = 1 solves 
Laplace's equation and Ui = 1 has differences equal to zero. The constant vector 
solves KU = except near the boundaries. When a neighbor is a boundary point 
where Ui is known, its value moves onto the right side of KU = F. Then that row of 
K is not zero sum. Otherwise K would be singular, if K * ones(n, 1) = zeros(?7,, 1). 

Using block matrix notation, we can create the 2D matrix K = K2D from the 
familiar by second difference matrix K. We number the nodes of the square a 



3.6. SOLVING LARGE LINEAR SYSTEMS 



193 



row at a time (this "natural numbering" is not necessarily best). Then the — I's for 
the neighbor above and the neighbor below are positions away from the main 
diagonal of K2d- The 2D matrix is block tri diagonal with tridiagonal blocks: 



K 



-1 
2 



-1 2 



K2D = 



K 



f 21 -I 
-I K + 2I 



-I 



Size N 
Time N 



-I K + 2I 

Elimination in this order: K2D has size n 



Bandwidth w = N, Space nw = N^, Time nw 



Ar2 
2 _ 



The matrix K2D has 4's down the main diagonal. Its bandwidth w = N is the 
distance from the diagonal to the nonzeros in — /. Many of the spaces in between are 
filled during elimination! Then the storage space required for the factors in K = LU 
is of order nw = N^. The time is proportional to nw"^ = N^, when n rows each 
contain w nonzeros, and w nonzeros below the pivot require elimination. 

Those counts are not impossibly large in many practical 2D problems (and we 
show how they can be reduced). The horrifying large counts come for K3D in three 
dimensions. Suppose the 3D grid is numbered by square cross-sections in the natural 
order 1, . . . , A^. Then K3D has blocks of order N"^ from those squares. Each square 
is numbered as above to produce blocks coming from K2D and I = I2D: 



K 



3D 



K2D + 21 -I 

-I K2D + 2I -I 



-I K2n + 2I 



Size n = 

Bandwidth w = N"^ 
Elimination space nw 
Elimination time k, nw 



ATS 



Now the main diagonal contains 6's, and "inside rows" have six —I's. Next to a point 
or edge or corner of the boundary cube, we lose one or two or three of those —I's. 

The good way to create K2D from K and / (A^ by A^) is to use the kron[A,B) 
command. This Kronecker product replaces each entry aij by the block aijB. To take 
second differences in all rows at the same time, and then all columns, use kron: 



K2D = kron{K, I) + kron(J, K) . 



(2) 



The identity matrix in two dimensions is I2D = kron (/,/). This adjusts to allow 
rectangles, with /'s of different sizes, and in three dimensions to allow boxes. For a 
cube we take second differences inside all planes and also in the ^-direction: 

KsD = kron(K2D, I) + kron(l2D, K) ■ 



Having set up these special matrices K2D and K3D, we have to say that there are 
special ways to work with them. The x, y, z directions are separable. The geometry 
(a box) is also separable. See Section 7.2 on Fast Poisson Solvers. Here the matrices 
K and K2D and K3D are serving as models of the type of matrices that we meet. 
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Minimum Degree Algorithm 

We now describe (a little roughly) a useful reordering of the nodes and the equations 
in K2dU = F. The ordering achieves minimum degree at each step — the number 
of nonzeros below the pivot row is minimized. This is essentially the algorithm 
used in MATLAB's command U = K\F, when K has been defined as a sparse matrix. 
We list some of the functions from the sparfun directory: 



speye (sparse identity /) nnz (number of nonzero entries) 

find (find indices of nonzeros) spy (visualize sparsity pattern) 

colamd and symamd (approximate minimum degree permutation of K) 



You can test and use the minimum degree algorithms without a careful analysis. The 
approximations are faster than the exact minimum degree permutations colmmd and 
symmmd. The speed (in two dimensions) and the roundoff errors are quite reasonable. 

In the Laplace examples, the minimum degree ordering of nodes is irregular com- 
pared to "a row at a time." The final bandwidth is probably not decreased. But the 
nonzero entries are postponed as long as possible! That is the key. 

The difference is shown in the arrow matrix of Figure 3. 18. On the left, minimum 
degree (one nonzero off the diagonal) leads to large bandwidth. But there is no fill-in. 
Elimination will only change its last row and column. The triangular factors L and 
U have all the same zeros as A. The space for storage stays at 3n, and elimination 
needs only n divisions and multiplications and subtractions. 



* 



Bandwidth 
6 and 3 

Fill-in 
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* 
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Figure 3.18: Arrow matrix: Minimum degree (no F) against minimum bandwidth. 



The second ordering reduces the bandwidth from 6 to 3. But when row 4 is 
reached as the pivot row, the entries indicated by F are filled in. That full lower 
quarter of A gives nonzeros to both factors L and U. You see that the whole 
"profile" of the matrix decides the fill-in, not just the bandwidth. 

The minimum degree algorithm chooses the {k + l)st pivot column, after k 
columns have been eliminated as usual below the diagonal, by the following rule: 

In the remaining matrix of size n — k, select the column with the fewest nonzeros. 
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The component of U corresponding to that column is renumbered k + 1. So is the 
node in the finite difference grid. Of course ehmination in that column will normally 
produce new nonzeros in the remaining columns! Some fill-in is unavoidable. So 
the algorithm must keep track of the new positions of nonzeros, and also the actual 
entries. It is the positions that decide the ordering of unknowns. Then the entries 
decide the numbers in L and U. 

Example Figure 3.19 shows a small example of the minimal degree ordering, for 
Laplace's 5-point scheme. The node connections produce nonzero entries (indicated by 
*) in K. The problem has six unknowns. K has two 3 by 3 tridiagonal blocks from 
horizontal links, and two 3 by 3 blocks with — / from vertical links. 

The degree of a node is the number of connections to other nodes. This is the 
number of nonzeros in that column of K. The corner nodes 1, 3, 4, 6 all have degree 
2. Nodes 2 and 5 have degree 3. A larger region has inside nodes of degree 4, which 
will not be eliminated first. The degrees change as elimination proceeds, because of 
fill-in. 

The first elimination step chooses row 1 as pivot row, because node 1 has minimum 
degree 2. (We had to break a tie! Any degree 2 node could come first, leading to different 
elimination orders.) The pivot is P, the other nonzeros in that row are boxed. When 
row 1 operates on rows 2 and 4, it changes six entries below it. In particular, the two 
fill-in entries marked byF change to nonzeros. This fill-in of the (2,4) and (4,2) entries 
corresponds to the dashed line connecting nodes 2 and 4 in the graph. 
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Figure 3.19: Minimum degree nodes 1 and 3. The pivots P are in rows 1 and 3; new 
edges 2-4 and 2-6 in the graph match the matrix entries F filled in by elimination. 
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Nodes that were connected to the eliminated node are now connected to each other. 
Elimination continues on the 5 by 5 matrix (and the graph with 5 nodes). Node 2 still has 
degree 3, so it is not eliminated next. If we break the tie by choosing node 3, elimination 
using the new pivot P will fill in the (2,6) and (6,2) positions. Node 2 becomes linked 
to node 6 because they were both linked to the eliminated node 3. 

The problem is reduced to 4 by 4, for the unknown U's at the remaining nodes 

2,4,5,6. Problem asks you to take the next step — choose a minimum degree node 

and reduce the system to 3 by 3. 

Storing the Nonzero Structure = Sparsity Pattern 

A large system KU = F needs a fast and economical storage of the node connections 
(which match the positions of nonzeros in K). The connections and nonzeros change 
as elimination proceeds. The list of edges and nonzero positions corresponds to the 
adjacency matrix" of the graph of nodes. The adjacency matrix has 1 or to indicate 
nonzero or zero in K. 

For each node i, we have a list adj(z) of the nodes connected to i. How to combine 
these into one master list NZ for the whole graph and the whole matrix K? A simple 
way is to store the lists a.d]{i) sequentially in NZ (the nonzeros for i = 1 up to i = n). 
An index array IND of pointers tells the starting position of the sublist adj(2) within 
the master list NZ. It is useful to give IND an (n + l)st entry to point to the final 
entry in NZ (or to the blank that follows, in Figure 3.20). MATLAB will store one 
more array (the same length nnz(K) as NZ) to give the actual nonzero entries. 

NZ |24|135|26|15|246|35| 

IND 1 3 6 8 10 13 15 

Node 1 2 3 4 5 6 

Figure 3.20: Master list NZ of nonzeros (neighbors in Figure 3.19). Positions in 
IND. 

The indices i are the "original numbering" of the nodes. If there is renumbering, 
the new ordering can be stored as a permutation PERM. Then PERM(2) = k when 
the new number i is assigned to the node with original number k. The text [GL] by 
George and Liu is the classic reference for this entire section on ordering of the nodes. 

Graph Separators 

Here is another good ordering, different from minimum degree. Graphs or meshes are 
often separated into disjoint pieces by a cut. The cut goes through a small number 
of nodes or meshpoints (a separator). It is a good idea to number the nodes in the 
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separator last. Elimination is relatively fast for the disjoint pieces P and Q. It only 
slows down at the end, for the (smaller) separator S. 

The three groups P,Q,Sof meshpoints have no direct connections between P and 
Q (they are both connected to the separator S). Numbered in that order, the "block 
arrow" stiffness matrix and its K = LU factorization look like this: 



K 



Kp 





Kps 







Kqs 


Ksp 


KsQ 


Ks 



Lp 


X 



Lq 
Y 
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Uo 



A 
B 
C 



(3) 

The zero blocks in K give zero blocks in L and U. The submatrix Kp comes first 
in elimination, to produce Lp and Up. Then come the factors LqUq of Kq, followed 
by the connections through the separator. The major cost is often that last step, the 
solution of a fairly dense system of the size of the separator. 




Arrow matrix 

(Figure 3.18) 



S 



6 
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Separator comes last 
(Figure 3.19) 
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Blocks P, Q 
Separator S 



Figure 3.21: A graph separator numbered last produces a block arrow matrix K. 



Figure 3.21 shows three examples, each with separators. The graph for a perfect 
arrow matrix has a one-point separator (very unusual). The 6-node rectangle has a 
two-node separator in the middle. Every A by A grid can be cut by an A-point 
separator (and A is much smaller than A^). If the meshpoints form a rectangle, the 
best cut is down the middle in the shorter direction. 

You could say that the numbering of P then Q then 5* is block minimum degree. 
But one cut with one separator will not come close to an optimal numbering. It 
is natural to extend the idea to a nested sequence of cuts. P and Q have their 
own separators at the next level. This nested dissection continues until it is not 
productive to cut further. It is a strategy of "divide and conquer." 

Figure 3.22 illustrates three levels of nested dissection on a 7 by 7 grid. The first 
cut is down the middle. Then two cuts go across and four cuts go down. Numbering 
the separators last within each stage, the matrix K of size 49 has arrows inside arrows 
inside arrows. The spy command will display the pattern of nonzeros. 

Separators and nested dissection show how numbering strategies are based on the 
graph of nodes and edges in the mesh. Those edges correspond to nonzeros in the 
matrix K. The nonzeros created by elimination (filled entries in L and U) correspond 
to paths in the graph. In practice, there has to be a balance between simplicity and 
optimality in the numbering — in scientific computing simplicity is a very good thing! 
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Figure 3.22: Three levels of separators. Still 
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A very reasonable compromise is the backslash command U = K\F that uses a 
nearly minimum degree ordering in Sparse MATLAB. 



Operation Counts (page K) 

Here are the complexity estimates for the 5-point Laplacian with A^^ or A^'^ nodes: 

Minimum Degree n = N"^ in 2D n = in 3D 

Space (nonzeros from fill-in) X X 

Time (flops for elimination) X X 

Nested Dissection 

Space (nonzeros from fill-in) X X 

Time (flops for elimination) X X 

In the last century, nested dissection lost out — it was slower on almost all applica- 
tions. Now larger problems are appearing and the asymptotics eventually give nested 
dissection an edge. Algorithms for cutting graphs can produce short cuts into nearly 
equal pieces. Of course a new idea for ordering could still win. 



Iterative versus Direct Methods 

This section is a guide to solution methods for problems Ax = b that are too large 
and expensive for ordinary elimination. We are thinking of sparse matrices A, when 
a multiplication Ax is relatively cheap. If A has at most p nonzeros in every row, 
then Ax needs at most pn multiplications. Typical applications are to large finite 
difference equations or finite element problems on unstructured grids. In the special 
case of a square grid for Laplace's equation, a Fast Poisson Solver (Section 7.2) is 
available. 

We turn away from elimination to iterative methods and Krylov subspaces. 

Pure iterative methods are easier to analyze, but the Krylov subspace methods are 
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more powerful. So the older iterations of Jacobi and Gauss-Seidel and overrelaxation 
are less favored in scientific computing, compared to conjugate gradients and GM- 
RES. When the growing Krylov subspaces reach the whole space R", these methods 
(in exact arithmetic) give the exact solution A~^b. But in reality we stop much ear- 
lier, long before n steps are complete. The conjugate gradient method {for positive 
definite A, and with a good preconditioner) has become truly important. 

The next ten pages will introduce you to numerical linear algebra. This has 
become a central part of scientific computing, with a clear goal: Find a fast stable 
algorithm that uses the special properties of the matrices. We meet matrices that 
are sparse or symmetric or triangular or orthogonal or tridiagonal or Hessenberg or 
Givens or Householder. Those matrices are at the core of so many computational 
problems. The algorithm doesn't need details of the entries (which come from the 
specific application). By using only their structure, numerical linear algebra offers 
major help. 

Overall, elimination with good numbering is the first choice until storage and CPU 
time become excessive. This high cost often arises first in three dimensions. At that 
point we turn to iterative methods, which require more expertise. You must choose 
the method and the preconditioner. The next pages aim to help the reader at this 
frontier of scientific computing. 

Pure Iterations 

We begin with old-style pure iteration (not obsolete). The letter K will be reserved 
for "Krylov" so we leave behind the notation KU = F. The linear system becomes 
Ax = b with a large sparse matrix A, not necessarily symmetric or positive definite: 

Linear system Ax = b Residual = b — Axk Preconditioner P ^ A 

The preconditioner P attempts to be "close to A" and at the same time much easier 
to work with. A diagonal P is one extreme (not very close). P = A is the other 
extreme (too close). Splitting the matrix A gives an equivalent form of Ax = b: 

Splitting Px = {P -A)x + b. (4) 

This suggests an iteration, in which every vector Xk leads to the next Xk+i'. 

Iteration Pxk+i = (P — A)xk + b . (5) 

Starting from any xq, the first step finds Xi from Pxi = {P — A)xo + b. The iteration 
continues to X2 with the same matrix P, so it often helps to know its triangular factors 
L and U. Sometimes P itself is triangular, or its factors L and U are approximations 
to the triangular factors of A. Two conditions on P make the iteration successful: 

1. The new Xk+i must be quickly computable. Equation (5) must be fast to solve. 

2. The errors = x — Xk must converge quickly to zero. 
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Subtract equation (5) from (4) to find the error equation. It connects to Ck+i'- 

Error Pck+i = (P — A)ek whicli means Ck+i = {I — P^^A)ek = Met ■ (6) 

Tlie riglit side h disappears in this error equation. Each step multiphes the error 
vector by M = / — P~^A. The speed of convergence of to x (and of to zero) 
depends entirely on M. The test for convergence is given by the eigenvalues of M: 

Convergence test Every eigenvalue of M must have |A(M)| < 1. 

The largest eigenvalue (in absolute value) is the spectral radius p{M) = max |A(M)|. 
Convergence requires p(M) < 1. The convergence rate is set by the largest eigen- 
value. For a large problem, we are happy with p{M) = .9 and even p{M) = .99. 

Suppose that the initial error cq happens to be an eigenvector of M. Then the 
next error is Ci = Mcq = Xcq. At every step the error is multiplied by A, so we must 
have |A| < 1. Normally cq is a combination of all the eigenvectors. When the iteration 
multiplies by M, each eigenvector is multiplied by its own eigenvalue. After k steps 
those multipliers are A'^. We have convergence if all |A| < 1. 

For preconditioner we first propose two simple choices: 

Jacobi iteration P = diagonal part of A 

Gauss-Seidel iteration P = lower triangular part of A 

Typical examples have spectral radius p(M) = 1 — cN~^. This comes closer and 
closer to 1 as the mesh is refined and the matrix grows. An improved preconditioner 
can give p(M) = 1 - cN-^/^. Then p is smaller and convergence is faster, as in 
"overrelaxation." But a different approach has given more flexibility in constructing 
a good P, from a quick incomplete LU factorization of the true matrix A: 

/ncomplete LU P = (approximation to L) (approximation to U) . 

The exact A = LU has fill-in, so zero entries in A become nonzero in L and U. The ap- 
proximate L and U could ignore this fill-in (fairly dangerous). Or P = i^approx^approx 
can keep only the fill-in entries F above a fixed threshold. The variety of options, 
and the fact that the computer can decide automatically which entries to keep, has 
made the ILU idea (incomplete LU) a very popular starting point. 

Example The —1,2,-1 matrix A = K provides an excellent example. We choose 
the preconditioner P = T, the same matrix with Tn = 1 instead of Kn = 2. The LU 
factors of T are perfect first differences, with diagonals of -|-1 and —1. (Remember that 
all pivots of T equal 1, while the pivots of K are 2/1, 3/2, 4/3, . . .) We can compute the 
right side of T^^Kx = T^^b with only 2A^ additions and no multiplications (just back 
substitution using L and U). Idea: This L and U are approximately correct for K. 

The matrix P~^A = T^^K on the left side is triangular. More than that, T is a rank 1 
change from K (the 1, 1 entry changes from 2 to 1). It follows that T~^K and K^^T 
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are rank 1 changes from the identity matrix /. A calculation in Problem shows that 

only the first column of I is changed, by the "linear vector" £ = {N, N — 1, . . . ,1): 

p-^A = T~^K = I + £ej and R-^T = I ~ {£eJ)/{N + 1) . (7) 

Here e]" = [ 1 ... ] so £ej has first column £. This example finds x = K~^b by 
a quick exact formula {K^^T)T~^h, needing only 2N additions for and N additions 
and multiplications for K~^T. In practice we wouldn't precondition this K (just solve). 

The usual purpose of preconditioning is to speed up convergence for iterative methods, 
and that depends on the eigenvalues of P^^A. Here the eigenvalues of T^^K are its 
diagonal entries A^+1, 1, . . . , 1. This example will illustrate a special property of conjugate 
gradients, that with only two different eigenvalues it reaches the true solution x in two 
steps. 

The iteration Pxk+i = {P — A)xk + b is too simple! It is choosing one particular 
vector in a "Krylov subspace." With relatively little work we can make a much better 
choice of Xk- Krylov projections are the state of the art in today's iterative methods. 

Krylov Subspaces 

Our original equation is Ax = b. The preconditioned equation is P"^Ax = P^^b. 
When we write P^^, we never intend that an inverse would be explicitly computed 
(except in our example). The ordinary iteration is a correction to Xk by the vector 

p-w. 

Pxk+i = {P - A)xk + b or Pxk+i = Pxk + Vk or Xk+i = Xk + P'^Vk . (8) 

Here Vk = b — Axk is the residual. It is the error in Ax = b, not the error Ck 
in X. The symbol P'^r^ represents the change from Xk to Xk+i, but that step is 
not computed by multiplying P~^ times r^. We might use incomplete LU , or a few 
steps of a "multigrid" iteration, or "domain decomposition." Or an entirely new 
preconditioner. 

In describing Krylov subspaces, I should work with P~^A. For simplicity I will 
only write A. I am assuming that P has been chosen and used, and the precondi- 
tioned equation P~^Ax = P~^b is given the notation Ax = b. The preconditioner is 
now P = I. Our new matrix A is probably better than the original matrix with that 
name. 

The Krylov subspace K^(A, b) contains all combinations of b, Ab, . . . , A''~^b. 

These are the vectors that we can compute quickly, multiplying by a sparse A. We 
look in this space K'^ for the approximation Xk to the true solution of Ax = b. Notice 
that the pure iteration Xk = {I — A)xk-i + b does produce a vector in K'^ when Xk-i 
is in K'^"^. The Krylov subspace methods make other choices of Xk- Here are four 
different approaches to choosing a good Xk in K'^ — this is the important decision: 
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1. The residual = 6 — Axk is orthogonal to K'^ (Conjugate Gradients, . . . ) 

2. The residual has minimum norm for Xk in K'^' (GMRES, MINRES, . . . ) 

3. Tfc is orthogonal to a different space like K.'^{A'^) (BiConjugate Gradients, . . . ) 

4. Cfc has minimum norm (SYMMLQ; for BiCGStab Xk is in A^K^'^A^); . . . ) 



In every case we hope to compute the new Xk quickly and stably from the earlier x's. 
If that recursion only involves Xk-i and Xk-2 (short recurrence) it is especially fast. 
We will see this happen for conjugate gradients and symmetric positive definite A. 
The BiCG method in 3 is a natural extension of short recurrences to unsymmetric 
A — but stability and other questions open the door to the whole range of methods. 

To compute x^ we need a basis for K'^. The best basis gi, . . . , is orthonormal. 
Each new comes from orthogonalizing t = Aq^-i to the basis vectors gi, . . . , qk-i 
that are already chosen. This is the Gram-Schmidt idea (called modified Gram- 
Schmidt when we subtract projections of t onto the g's one at a time, for numerical 
stability). The iteration to compute the orthonormal g's is known as Arnoldi's 
method: 



1 qi = b/\\b\\2; % Normalize to ||gi|| = 1 
for j = 1, . . . , A; — 1 

2 t = Aqj; % t is in the Krylov space K^+^{A, b) 
for z = 1, . . . , j 

3 hij = qjt] % hijqi = projection of t onto qi 

4 t = t — hijqi, % Subtract component of t along qi 
end; 

5 ^i+i.i = ll^lh; % t is now orthogonal to gi, . . . , g^ 

6 gj+i = t/hj+ij] % Normalize t to ||gj+i|| = 1 

end % gi, . . . , gfc are orthonormal in K'^ 



Put the column vectors gi, . . . , g^ into an n by /c matrix Qk- Multiplying rows of 
Qj by columns of Qk produces all the inner products g^^gj, which are the O's and I's 
in the identity matrix. The orthonormal property means that QlQk = Ik. 

Arnoldi constructs each g^+i from Aqj by subtracting projections hijqi. If we 
express the steps up to j = A; — 1 in matrix notation, they become AQk-i = QkHk^k-i'- 



Arnoldi 

AQk-i 



Aqi 



Mk- 



Qk 



n hj k — 1 



nhj k 



hii hi2 ■ hi^k-i 

^21 ^22 ■ ^2,A:-l 

Q h2z • 

• hk,k-i 

khy k — 1 



(9) 



That matrix Hk^k-i is "upper Hessenberg" because it has only one nonzero diagonal 
below the main diagonal. We check that the first column of this matrix equation 
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(multiplying by columns!) produces q2. 

Aqi = hiiQi + h2iq2 or q2 = r . (10) 

'^21 

That subtraction is Step 4 in Arnoldi's algorithm. Division by /121 is Step 6. 

Unless more of the hij are zero, the cost is increasing at every iteration. We have 
k dot products to compute at step 3 and 5, and k vector updates in steps 4 and 6. A 
short recurrence means that most of these hij are zero. That happens when A = A"^. 

The matrix H is tridiagonal when A is symmetric. This fact is the founda- 
tion of conjugate gradients. For a matrix proof, multiply equation (9) by Q^-i- The 
right side becomes H without its last row, because {Ql_iQk)Hk,k~i = [I 0] Hk k_i. 
The left side Qj._]^AQk^i is always symmetric when A is symmetric. So that H matrix 
has to be symmetric, which makes it tridiagonal. There are only three nonzeros in 
the rows and columns of if, and Gram-Schmidt to find qk+i only involves qk and qt-i'. 

Arnoldi when A = Aqt = hk+i,k qk+i + hk,k qk + hk-i,k qk-i ■ (H) 

This is the Lanczos iteration. Each new qk+i = {Aqk — hk,kqk — hk^i,kqk-i) / hk+i,k 
involves one multiplication Aq^, two dot products for new /I's, and two vector updates. 

The QR Method for Eigenvalues 

Allow me an important comment on the eigenvalue problem Ax = Ax. We have seen 
that Hk_i = Q'^_^AQk-i is tridiagonal if A = A^. When k — 1 reaches n and Qn is 
square, the matrix H = Q^AQ^ = Q~^AQn has the same eigenvalues as A: 

Same A Hy = Q~^AQny = \y gives Ax = Xx with x = Qny • (12) 

It is much easier to find the eigenvalues A for a tridiagonal H than the for original A. 

The famous "QR method" for the eigenvalue problem starts with Ti = H, factors 
it into Ti = QiRi (this is Gram-Schmidt on the short columns of Ti), and reverses 
order to produce T2 = RiQi. The matrix T2 is again tridiagonal, and its off-diagonal 
entries are normally smaller than for Ti. The next step is Gram-Schmidt on T2, 
orthogonalizing its columns in Q2 by the combinations in the upper triangular R2: 

QR Method Factor T2 into Q2R2 ■ Reverse order to T3 = R2Q2 = ^^2<52 (13) 

By the reasoning in (12), any Q^^TQ has the same eigenvalues as T. So the matrices 
T2,T3, . . . all have the same eigenvalues as, Ti = H and A. (These square Qk from 
Gram-Schmidt are entirely different from the rectangular Qk in Arnoldi.) We can 
even shift T before Gram-Schmidt, and we should, provided we remember to shift 
back: 

Shifted QR Factor Tk — Skl = QkRk ■ Reverse to Tk+i = RkQk + bkl ■ (14) 

When the shift Sk is chosen to be the ?T,,n entry of T^, the last off-diagonal entry of 
T/;_|_i becomes very small. The ra, n entry of T^+i moves close to an eigenvalue. Shifted 
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QR is one of the great algorithms of numerical linear algebra. It solves moderate-size 
eigenvalue problems with great efficiency. This is the core of MATLAB's e\g{A). 

For a large symmetric matrix, we often stop the Arnoldi-Lanczos iteration at a 
tridiagonal Hk with k < n. The full n-step process to reach Hn is too expensive, and 
often we don't need all n eigenvalues. So we compute (by the same QR method) the k 
eigenvalues of instead of the n eigenvalues of Hn. These computed Ai^, A2A;, . . . , Xkk 
can provide good approximations to the first k eigenvalues of A. And we have an 
excellent start on the eigenvalue problem for Hk+i, if we decide to take a further step. 

This Lanczos method will find, approximately and iteratively and quickly, the 
leading eigenvalues of a large symmetric matrix. 

The Conjugate Gradient Method 

We return to iterative methods for Ax = b. The Arnoldi algorithm produced or- 
thonormal basis vectors qi,q2, . . . for the growing Krylov subspaces K^, K^, . . .. Now 
we select vectors Xi,X2,... in those subspaces that approach the exact solution to 
Ax = b. We concentrate on the conjugate gradient method for symmetric positive 
definite A. 

The rule for Xk in conjugate gradients is that the residual = b — Ax^ should 
be orthogonal to all vectors in K'^. Since will be in K'^"'"^, it must be a multiple 
of Arnoldi's next vector gfc+i! Each residual is therefore orthogonal to all previous 
residuals (which are multiples of the previous g's): 

Orthogonal residuals rfr^ = for i < k . (15) 

The difference between and qk+i is that the g's are normalized, as in qi = b/\\b\\. 

Since r^.i is a multiple of q^, the difference r^— r^.i is orthogonal to each subspace 
K* with i < k. Certainly Xi — Xi-i lies in that K*. So Ar is orthogonal to earlier 
Ax's: 

{xi - Xi_i)'^(rfc - rfc_i) =0 for z < . (16) 
These differences Ax and Ar are directly connected, because the fe's cancel in Ar: 

rfc - rfc_i = (6 - Axk) - (6 - Axk-i) = -A{xk - Xk-i) . (17) 

Substituting (17) into (16), the updates in the x's are "A-orthogonal" or conjugate: 

Conjugate updates Ax (xj — Xi-i)'^A(xk — x^-i) = for i < k. (18) 

Now we have all the requirements. Each conjugate gradient step will find a new 
"search direction" dk for the update Xk — Xk-i. From Xk-i it will move the right 
distance akdk to Xk. Using (17) it will compute the new r^. The constants Pk in 
the search direction and ak in the update will be determined by (15) and (16) for 
i = k — 1. For symmetric A the orthogonality in (15) and (16) will be automatic for 
z < A; — 1, as in Arnoldi. We have a "short recurrence" for the new x^ and r^. 
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Here is one cycle of the algorithm, starting from xq = and tq = b and /3i = 0. It 
involves only two new dot products and one matrix-vector multiplication Ad: 



Conjugate 1 Pk = ^k-i^k-i/i^k-2''^k-2 % Improvement this step 

Gradient 2 dk = r^-i + (3kdk-i % Next search direction 

Method 3 = r'^_^rk-i/djAdk % Step length to next Xk 

4 Xk = Xk-i + akdk % Approximate solution 

5 rk = Tk-i — akA.dk % New residual from (17) 



The formulas 1 and 3 for j3k and ak are explained briefly below — and fully by Trefethen- 
Bau ( ) and Shewchuk ( ) and many other good references. 

Different Viewpoints on Conjugate Gradients 

I want to describe the (same!) conjugate gradient method in two different ways: 

1. It solves a tridiagonal system Hy = f recursively 

2. It minimizes the energy ^x"^Ax — x'^b recursively. 

How does Ax = b change to the tridiagonal Hy = /? That uses Arnoldi's or- 
thonormal columns gi, . . . , g„ in Q, with Q^Q = I and Q^AQ = H: 

Ax = b is {Q^AQ){Q^x) = Q^b which is //y = / = (||6||, 0, . . . , 0) . (19) 

Since qi is the first component of / = Q^b is qjb = \\b\\ and the other com- 

ponents are qjb = 0. The conjugate gradient method is implicitly computing this 
symmetric tridiagonal H and updating the solution y at each step. Here is the third 
step: 
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hl2 
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h23 




y3 











hz2 


^33 












This is the equation Ax = b projected by Q3 onto the third Krylov subspace K^. 

These /I's never appear in conjugate gradients. We don't want to do Arnoldi too! 
It is the LDL^ factors of H that CG is somehow computing — two new numbers at 
each step. Those give a fast update from yj^i to yj. The corresponding Xj = Qjyj 
from conjugate gradients approaches the exact solution Xn = QnVn which is x = A~^b. 

If we can see conjugate gradients also as an energy minimizing algorithm, we can 
extend it to nonlinear problems and use it in optimization. For our linear equation 
Ax = b, the energy is E{x) = \x^ Ax — x^b. Minimizing E{x) is the same as solving 
Ax = b, when A is positive definite (the main point of Section 1. ). The CG 
iteration minimizes E{x) on the growing Krylov subspaces. On the first 
subspace K^, the line where x is ab = adi, this minimization produces the right 



206 



CHAPTER 3. BOUNDARY VALUE PROBLEMS 



value for ai: 



E{ah) 



1 



a^l)^ Ah — ab^b is minimized at ai 



b^b 



(21) 



That ai is the constant chosen in step 3 of the first conjugate gradient cycle. 



The gradient of E{x) 



_1 A , 

2 ■X' jj.iAj 



X 



b is exactly Ax — b. The steepest descent 



direction at xi is along the negative gradient, which is ri! This sounds like the perfect 
direction d2 for the next move. But the great difficulty with steepest descent is that 
this Ti can be too close to the first direction. Little progress that way. So we add 
the right multiple P2di, in order to make d2 = ri + P2di A-orthogonal to the first 
direction di. 

Then we move in this conjugate direction d2 to X2 = Xi + 02^2- This explains the 
name conjugate gradients, rather than the pure gradients of steepest descent. Every 



X 



cycle of CG chooses aj to minimize E{x) in the new search direction x 
The last cycle (if we go that far) gives the overall minimizer Xr. 

Example 



x^i + 



adj. 



A"^b. 



Ax 
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'211" 
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From Xo = and f3i = and = di = b the first cycle gives ai = \ and Xi 
(2, 0, 0). The new residual \s ri = b — Axi = (0, —2, —2). Then the second cycle yields 



^^-16' 



do 



2 
-2 
-2 



"^=16' 



X2 



3 
-1 
-1 



A-^b\ 



The correct solution is reached in two steps, where normally it will take n = 3 steps. The 
reason is that this particular A has only two distinct eigenvalues 4 and 1. In that case 
A~^b is a combination of 6 and Ab, and this best combination X2 is found at cycle 2. The 
residual r2 is zero and the cycles stop early — very unusual. 

Energy minimization leads in [ ] to an estimate of the convergence rate for the 



error e 



X 



X, 



in conjugate gradients, using the A- norm \\e\\A = V e^Ae: 



Error estimate \\x — Xj\\A < 2 ( ^^^^^ — V-^mm \ \]^^ _ ^^\^\^^ _ (^22) 

This is the best-known error estimate, although it doesn't account for any clustering of 

the eigenvalues of A. It involves only the condition number Amax/ Amin- Problem 

gives the "optimal" error estimate but it is not so easy to compute. That optimal 
estimate needs all the eigenvalues of A, while (22) uses only the extreme eigenvalues 
Amax(^) and Amin(v4) — which in practice we can bound above and below. 
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Minimum Residual Methods 

When A is not symmetric positive definite, conjugate gradient is not guaranteed to 
solve Ax = b. Most likely it won't. We will follow van der Vorst [ ] in briefly 
describing the minimum norm residual approach, leading to MINRES and GMRES. 

These methods choose Xj in the Krylov subspace K-' so that — AxjH is minimal. 
First we compute the orthonormal Arnoldi vectors qi, . . . ,qj. They go in the columns 
of Qj, so QjQj = I. As in (19) we set Xj = QjU, to express the solution as a 
combination of those g's. Then the norm of the residual Vj using (9) is 

II 6 - AxjW = \\b- AQjyW = \\b - Qj+iHj+ijy\\ . (23) 

These vectors are all in the Krylov space K-'^^, where rJ{Qj^iQj_^^rj) = rjrj. This 
says that the norm is not changed when we multiply by Qj^i- Our problem becomes: 

Choose y to minimize ||rj|| = WQj^ib — Hj^ijy\\ = \\f — Hy\\ . (24) 

This is an ordinary least squares problem for the equation Hy = f with only j + 1 
equations and j unknowns. The right side / = Qj^ib is (||ro||,0, . . . , 0) as in (19). 
The matrix H = -ffj+ij is Hessenberg as in (9), with one nonzero diagonal below the 
main diagonal. We face a completely typical problem of numerical linear algebra: 
Use the special properties of H and f to find a fast algorithm that computes y. The 
two favorite algorithms for this least squares problem are closely related: 

MINRES A is symmetric (probably indefinite, or we use CG) and H is tridiagonal 
GMRES A is not symmetric and the upper triangular part of H can be full 

In both cases we want to clear out that nonzero diagonal below the main diagonal of 
H. The natural way to do that, one nonzero entry at a time, is by ^^Givens rotations." 
These plane rotations are so useful and simple (the essential part is only 2 by 2) that 
we complete this section by explaining them. 

Givens Rotations 

The direct approach to the least squares solution of Hy = f constructs the normal 
equations H"^Hy = H^f. That was the central idea in Chapter 1, but you see what 
we lose. If H is Hessenberg, with many good zeros, H^H is full. Those zeros in H 
should simplify and shorten the computations, so we don't want the normal equations. 

The other approach to least squares is by Gram-Schmidt. We factor H into 
orthogonal times upper triangular. Since the letter Q is already used, the or- 
thogonal matrix will be called G (after Givens). The upper triangular matrix is G^^H . 
The 3 by 2 case shows how a plane rotation G21 can clear out the subdiagonal entry 
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l'21- 



G21H 



COS 6 sin 6 
sin 6 cos 9 
1 



hii 


hi2 




* 


* 




h22 







* 





hz2 







* 



(25) 



That bold zero entry requires hu sin 9 = /121 cos 9, which determines 9. A second 
rotation G'^2^ in the 2-3 plane, will zero out the 3,2 entry. Then Gjg^Gg/i/ is a 
square upper triangular matrix U above a row of zeros! 

The Givens orthogonal matrix is G = G21G32 but there is no reason to do this mul- 
tiplication. We use each Gij as it is constructed, to simplify the least squares problem. 
Rotations (and all orthogonal matrices) leave the lengths of vectors unchanged: 



\\Hy - /II = \\G^2'G2i'Hy - G^iG^.'fW 



' u ' 




' F ' 





y - 


e 



(26) 



This length is what MINRES and GMRES minimize. The row of zeros below U 
means that the last entry e is the error — we can't reduce it. But we get all the other 
entries exactly right by solving the j by j system Uy = F (here j = 2). This gives 
the best least squares solution y. Going back to the original problem of minimizing 
Ik II — ~ ^^jll' ^he best Xj in the Krylov space K-' is Qjy. 

For non-symmetric A (GMRES rather than MINRES) we don't have a short 
recurrence. The upper triangle in H can be full, and step j becomes expensive and 
possibly inaccurate as j increases. So we may change "full GMRES" to GMRES (m), 
which restarts the algorithm every m steps. It is not so easy to choose a good m. 



Problem Set 3.6 



1 Create K2D for a 4 by 4 square grid with A^^ = 3^ interior mesh points (so 
n = 9). Print out its factors K = LU (or its Cholesky factor C = chol(K) for 
the symmetrized form K = C"'"C). How many zeros in these triangular factors? 
Also print out inv(K) to see that it is full. 

2 As increases, what parts of the LU factors of K2D are filled in? 

3 Can you answer the same question for K3D? In each case we really want an 
estimate cN^ of the number of nonzeros (the most important number is p). 

4 Use the tic; toe clocking command to compare the solution time for K2DX = 
random f in ordinary MATLAB and sparse MATLAB (where K2D is defined as 
a sparse matrix). Above what value of does the sparse routine K\f win? 

5 Compare ordinary vs. sparse solution times in the three-dimensional Kg^x = 
random f . At which does the sparse K\f begin to win? 

6 Incomplete LU 

7 Conjugate gradients 
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8 Draw the next step after Figure 3.19 when the matrix has become 4 by 4 and 
the graph has nodes 2-4-5-6. Which have minimum degree? Is there more 
fill-in? 

9 Redraw the right side of Figure 3.19 if row number 2 is chosen as the second 
pivot row. Node 2 does not have minimum degree. Indicate new edges in the 
5-node graph and new nonzeros F in the matrix. 

10 To show that T~^K = / + lef in (7), with e^'" = [ 1 ... ], we can start from 
K = T + eiej. Then T~^K = / + {T~'^ei)el and we verify that ei = Tt. 
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= ei . 



Second differences of a linear vector i are zero. Now multiply T = I + iej 
times / — {ieJ)/{N + 1) to establish the inverse matrix K~^T in (7). 

11 Arnoldi expresses each Aq^ as h^+i^k^k+i + h^^klk + ■ — I- Multiply by qj 

to find hi^k = Qi ^Qk- U ^ ^■^ symmetric you can write this as {Aqi)'^qk. Explain 

why {Aqi)'^qk = for z < A; - 1 by expanding Aqi into /li+i^^gi+i H h h^iqi. 

We have a short recurrence if A = A'^ (only hk+i,k and hk^k and hk-i,k are 
nonzero). 



